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Cationic lipid membranes ai'e known to form compact complexes with DNA and to be effective as gene delivery 
agents both in vitro and in vivo. Here we employ molecular dynamics simulations for a detailed atomistic study 
of lipid bilayers consisting of a mixture of cationic dimyristoyltrimethylammonium propane (DMTAP) and zwit- 
terionic dimyristoylphosphatidylcholine (DMPC). Our main objective is to examine how the composition of the 
bilayers affects their structural and electrostatic properties in the liquid-crystalline phase. By varying the mole 
fraction of DMTAP, we have found that the area per lipid has a pronounced non-monotonic dependence on the 
DMTAP concentration, with a minimum around the point of equimolar mixture. We show that this behavior has 
an electrostatic origin and is driven by the interplay between positively charged TAP headgroups and the zwitte- 
rionic PC heads. This inteiplay leads to considerable re-orientation of PC headgroups for an increasing DMTAP 
concentration, and gives rise to major changes in the electrostatic properties of the lipid bilayer, including a sig- 
nificant increase of total dipole potential across the bilayer and prominent changes in the ordering of water in the 
vicinity of the membrane. Moreover, chloride counter-ions are bound mostly to PC nitrogens implying stronger 
screening of PC heads by CI ions compared to TAP head groups. The implications of these findings are briefly 
discussed. 
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I. INTRODUCTION 

Gene therapy based on the introduction of genetic mate- 
rial into cells is one of the most promising biomedical ap- 
proaches t o treat human diseases jde Gennesl Il999t iLangeii 
fl998: .Monkkonen and Urttil Il998h . The majority of the de- 
livery vectors proposed are of viral nature. The viral vec- 
tors have been demonstrated to be very efficient but their 
use is restricted by accompanied toxicity ( Anderson. 1998) . 
This has stimulated a search for non-viral delivery sys- 
tems which should be characterized by greater safety and 
ease of manufacturing (Langer, 1998). Numerous exam- 
ples of non-viral delivery vectors include cationic liposomes, 
cationic polymers (such as polyamidoamine dendrimers, 
polyethylenimine, and spermine), and bl ock copolymers 
(p^tafieva et al., 1996; Fischer etaL, 1999; Gao and Huang, 
[1996; Ku kowskaLa tallo et al ., .1996|; Monkkonen and Urttit 
[1998; Pita rdetall ll997: Sm edt et alil20QCl) . 

In the light of the above, it is surprising how lit- 
tle attention has been devoted to computational studies 
of me mbranes containing cationic lipids. Bandyopadhyay 
et al. ('Bandvopadh vav et all '1999') performed an atomistic 
molecular dynamics (MD) study of a mixture of dimyris- 
toylphosphatidylcholine (DMPC) and dimyristoyltrimethy- 
lammonium propane (DMTAP) in the presence of a short 
DNA fragment. Apart from the very elegant piece of work 
above there are, to the best of our knowledge, no pub- 
lished atomistic computational studies of systems containing 



cationic lipids - this is very much in contrast to the great 
number of computational studies o f various neutral and an 
ionic phospholipid bilayer systems ("Fel leii 120001 ISaiz et al 



Tobias , 



2002[1Saiz and Klein. 2002; Tieleman et all 119971 
20011) . Another related example is the recent mo lecular dy- 
namics study of Bockmann et al. (iBockmann et alll2003) who 
showed the importance of monovalent ions on the properties 
and organization of lipid membranes - ions are always present 
in cationic lipid systems. The above examples demonstrate 
that detailed molecular dynamics studies can provide valuable 
insight into the atomistic organization of systems containing 
cationic lipids and yield useful information for experimen- 
talists about the underlying mechanisms on the atomic and 
molecular levels. 

In the present work, our objective is to gain insight into 
the structural and electrostatic properties of cationic lipid bi- 
layers through atomic classical molecular dynamics simula- 
tions. We concentrate on a bilayer mixture composed of two 
kinds of lipids: Neutral (zwitterionic) DMPC and cationic 
DMTAP (see Fig.[2for their chemical structures). Since DM- 
TAP is positively charged under physiological conditions, we 
have neutralized its positive charges by chloride counter-ions. 
From the computational point of view, this choice for a model 
system is motivated by the fact that DMPC and DMTAP have 
the same non-polar hydrocarbon chains and differ only by 
their headgroups. On the practical side, DMPC/DMTAP bi- 
nary lipid mixtures have been widely studied i n the presence 
of DNA by various experimental techniques dArtzner et all 
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FIG. 1 Chemical structures of the two lipids considered in 
the present work: a zwitterionic dimyristoylphosphatidylcholine 
(DMPC) and a cationic dimyristoyltrimethylammonium propane 
(DMTAP). 



Il998t iPohle et all |2000'; 'Zantl et al." 'l999a|d) and also by 
computational methods (.Bandvopadhvav et aUll99^ . 

We focus mostly on how the composition of the cationic 
lipid bilayer affects the structural and electrostatic properties 
of these lipid bilayer systems. To this end, we consider mix- 
tures of DMPC and DMTAP with various different mole frac- 
tions of the cationic DMTAP component under conditions cor- 
responding to the liquid-crystalline phase. We find that DM- 
TAP plays a prominent role leading to considerable changes 
compared to the pure DMPC bilayer As discussed in the 
present study, this is characterized by the strong interplay be- 
tween electrostatics and structural changes in the vicinity of 
the membrane-water interface. In particular, we find that DM- 
TAP gives rise to a non-monotonic dependence of the aver- 
age area per lipid on the DMTAP concentration, substantial 
changes in the electrostatic profile of the membrane, and sig- 
nificant re-orientation of P - N dipole vectors in DMPC head- 
groups. The spatial re-arrangement of P-N dipoles is par- 
ticularly interesting as it likely plays a significant role in the 
stability of DNA-membrane complexes. 



II. SYSTEM 

A. Model and simulation details 

We have performed atomistic simulations of fully hydrated 
lipid bilayers consisting of a mixture of cationic DMTAP and 
zwitterionic DMPC lipids. In all simulations, the total number 
of lipid molecules was fixed to 128, and the number of water 
molecules ranged from 3527 (pure DMTAP) to 3655 (pure 
DMPC). 

Force-field parameters for the lipids were taken from the re- 
cent united atom force-field (Ber ger et allll997|). This force- 
field has been previously validated ("Lin dahl and EdholmL 
I2OOO; Tieleman and Berendsen, 1996) and is essentially 
based on the GROMOS forcefield for lipid head groups, the 
Ryck a ert-Bel lman's potential jRvckaert and Bellemansl 
1 19751 Il978h for hydrocarbon chains and the OPLS 
parameters dJorgensen and Tirado-Rive^ Il988h for the 
Lennard-Jones interactions between united CH„ groups 



of acyl chains re-parameterized for long hydrocarbon 
chains to reproduce experimentally observed values 
of volume per hpid (|Vagle andWieneE Il988t) . The 
parameters for this force field are available on-line at 
http://moose.bio.ucalgary.ca/Downloads/files/lipid.itp Water 
was modeled using the SPC water model ( Berendsen et alj 
il981i) . The unit positive charge carried by each DMTAP 
molecule is compensated by the introduction of the cor- 
responding number of explicit Cl^ counter-ions. While 
being aware of the effect s of different models for chlo- 
ride dPatra and K arttuneru |200^, we decided to use the 
default set of chlor ide parameters suppli e d within the 
Grom acs force field (iBerendsen et all Il995t iLindahl et all 
I2OQII). 



Fol lowing the original parameterization ( Berge r et all 
11997b . the Lennard-Jones interactions were cut off at 1 nm 
without shift or switch function. Since long-range elec- 
trostatic interactions are essential in the present study, and 
since truncation of these interactions has been shown to 
lead to artifacts in simulations of phospholipids bilay- 
ers (Patraetal., 2003 ), we em ploy the particle-mesh Ewald 
(PME) method (.Darden et all Il993i) . The long-range contri- 
bution to the electrostatics is updated every 10-th time step. 

The simulations were performed in the NpT ensemble. 
The temperature was ke pt constant using a Berendsen ther- 
mostat (Berends en et all ll984) with a coupling time constant 
of 0.1 ps. Lipid molecules and water (including counter-ions) 
were separately coupled to a h eat bath. Pressure was con- 
trolled by a Berendsen barostat dBerendsen et all 1 1984ft with 
a coupling time constant of l.Ops. Pressure coupling was ap- 
plied semi-isotropically: The extension of the simulation box 
in the z direction (i. e., in the direction of the bilayer normal) 
and the cross-sectional area of the box in the x-y plane were 
allowed to vary independently of each other. 

We considered 1 1 DMPC / DMTAP mixtures ranging from 
pure DMPC to pure DMTAP. The molar fractions of the 
cationic DMTAP, xtap, were taken to be 0.0, 0.06, 0. 16, 0.25, 
0.31, 0.39, 0.50, 0.63, 0.75, 0.89, and 1.0. 

The main transition temperature of a pure DMPC bilayer is 
r„ = 24 °C dCevc and Mai'sh. .1987.). For DMPC / DMTAP 
binary mixtures it has been found Jzantl et allll999lJ) that the 
main transition temperature changes non-monotonic ally with 
the mole fraction of DMTAP, demonstrating a maximum of 
approximately 37 °C at xtap — 45. AU our simulations 
were done at a temperature of 50 °C, such that the bilayers are 
in the liquid-crystalline phase. 

All bond lengths of the lipid molecules were con- 
strained t o their equ i librium values using the LINCS al- 
gorithm iHess et al 1 Il997h whereas the SETTLE algo- 
rithm ( Mivamo to and Kollmani Il992h was used for water 
molecules. The time step in all simulations was set to 
2fs. All simulatio ns were perform ed using the GRO- 
MACS package (Ber endsen et allll995l:lLindahl et allEoOll) . 
The combined simulated time of all simulations amounts to 
250 nanoseconds. Each simulation was run in parallel over 
4 processors on an IBM eServer Cluster 1600 system. In total, 
the simulations took about 20,000 hours of CPU time. 
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B. Simulation setup 

Mixtures of DMPC and DMTAP were prepared and equili- 
brated in several steps as follows: 

1. We used the equilibrated DPP C bilayer structure of 
Patra et al. llPatra et all l2003h as our initial config- 
uration. The structure is available electronically at 
^http://www.softsimu.org/downloads.shtml_ 

2. DPPC and DMPC molecules differ only by length of 
their tail. Thus, we created a DMPC bilayer by short- 
ening the DPPC acyl chains by two hydrocarbons. This 
procedure does not disturb the acyl chain region or the 
water-hpid interface. 

3. The next step was to compress the DMPC bilayer in 
order to eliminate the gap between leaflets created in 
the previous step. For this purpose a pre-equilibration 
run for 1 ns in the NpT ensemble was performed. Af- 
ter this, the gap between leaflets had disappeared, and 
the obtained DMPC bilayer structure was used as ini- 
tial configuration for all DMPC / DMTAP mixtures de- 
scribed in the following. 

4. The chemical structures of DMPC and DMTAP dif- 
fer only by their head groups, see Fig. ^ With that 
in mind, the following procedure was used to prepare 
mixtures. For each DMTAP concentration the corre- 
sponding number of randomly chosen PC headgroups 
in a pu re DMPC bilayer were conv erted to TAP head 
groups (f Bandvopadhvav et alill999l) . To neutralize the 
unit charges in DMTAP headgroups, randomly chosen 
water molecules were replaced by chloride ions while 
ensuring a minimum separation of 0.5 nm between ions. 
To retain symmetry between the two leaflets, each of 
them contains the same number of cationic DMTAP. 

5. Since TAP headgroups occupy a smaller volume than 
fliose of DMPCs, a short lOps run in the NVT en- 
semble was performed to let the water molecules adjust 
at the lipid-water interface. This step completes pre- 
equilibration. 

6. The actual equilibration runs were performed in the 
NpT ensemble. The needed equilibration times prior 
to actual production runs ranged from 10 ns to 20 ns de- 
pending on the DMTAP mole fraction. We concluded 
that equilibration was completed when the average area 
per lipid had become stable and fluctuated around its 
mean with a standard deviation not exceeding the stan- 
dard deviation for a pure DMPC bilayer 

After equihbration, for each DMPC / DMTAP mixture, we 
performed a production run of 10 ns in the NpT ensemble to 
collect the data for analysis. 




Time [ns] 

FIG. 2 Time evolution of the area per lipid, A(t), for differ- 
ent mixtures of DMPC and DMTAP. The low concentration results 
(Xtap < 0.31) are shown on the top, and the high concentration 
regime (xtap > 0.63) is illustrated at the bottom. For clarity's 
sake, the results at intermediate concentrations are not shown here. 



III. RESULTS AND DISCUSSION 
A. Area per lipid 

The average area per lipid, {A), is one of the 
most fundamental characteristics of lipid bilay- 
ers (iNagle and Tristram- Naglei 120001) . While being one 
of the rather few structural quantities that can be measured 
accurately from model membranes via experiments, its also 
plays a major role in a number of quantities, including the 
ordering of acyl chains and the dynamics of lipids in a bilayer 
Further, from computational point of view, it is highly useful 
as a means of monitoring the equilibration process. 

Due to the lack of experimental data for the average area 
per lipid in binary DMPC / DMTAP bilayer mixtures, or in a 

J iure DMTAP bilay er [apart from the low-temperature phase 
Lewis et all l200ll) 1. reproduction of the experimental data 
available for pure DMPC membranes is essential to validate 
our approach. To this end, let us first consider the tempo- 
ral behavior of the area per lipid, A{t), presented in Fig. |2] 
It shows that the obtained average area per lipid for a pure 
DMPC system has a value oi {A) = 0.656 ± 0.008 nm^, in 
very good agreement with the experimentally observed value 
of 0.654 nm^ (.Petrache et all [2000 ). thereby vahdating our 
mode l in this respect. A s for the pure DMTAP bilayer, Lewis 
et al. llLewis et al .'.'200l') have been able to extract the area per 
lipid in the low- temperature phase, finding (A) = 0.40 nm^ at 
25° C. Studies of {A ) above the main transition temperature 
are lacking, however. 

We find that the average area per lipid shows a non- 
monotonic dependence on DMTAP mole fraction xtap, with 
a pronounced minimum roughly at xtap — 0.5, see Fig.|3l 
This behavior is not trivial, as modest amounts of the cationic 
DMTAP lead to a compression of the bilayer, while high con- 
centrations lead to a major expansion of the membrane. More 
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FIG. 3 Average area per lipid (A) as a function of the DMTAP 
mole fraction, xtap- 



0.72 

0.70 

0.68 

" 0.66 

I 0.64 

~^0.62 
A 

0.60 
0.58 
0.56 
0.54 



1 1 1 1 1 1 


1 1 1 1 1 1 1 1 1 1 1 


1 1 1 


- 


• DM PC 

DMTAP 


i 


5 ~- 






: i 












' 





III" 



0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 



FIG. 4 Average area per lipid {A)v based on the Voronoi analysis 
in two dimensions. The results for DMPC and DMTAP are shown 
separately as a function of xtap • 



specifically we find that for < xtap ^ 0.8, the average 
area per lipid is smaller than the corresponding counterpart 
for any of the pure lipid systems. Such a behavior cannot be 
explained by steric interaction alone but most likely is rather 
of electrostatic origin. 

While no published experimental data exists (to the best of 
our knowledge) for the area per lipid of DMPC /DMTAP bi- 
layer mixtures, similar behavior has been observed in related 
systems. 

For example, Zantl et al. (IZantl et all Il999bl) consid- 
ered DMPC /DMTAP monolayers using Langmuir-type film 
balance and found that, for pressures corresponding to the 
liquid-crystalline phase, the headgroup area decreased mono- 
tonically for small xtap, then had a minimum at about the 
equimolar ratio (xtap ~ 0.5), and increased for larger DM- 
TAP mole fractions. This is in accord with our observations. 

Another related work concerns Langmuir balance studies 
of mixed monolayers of zwitterionic palmitoyloleoylphos- 
phatidylcholine (POPC) and cationic 2,3-dimethoxy- 
l,4-bis(A^-hexadecy l-A^;A^-diniethyl-a mmonium)butane 
dibromide (SS-1) dSailv et all l200ll) . Although SS-1 is 
dicationic, one can qualitatively compare this system to 
our DMPC / DMTAP m i xture. For the POPC /SS-1 system 
Saily et al. dSailv et all 1200 ll) found that the average area 
per lipid has a non-monotonic behavior with a minimum at 
XSS-i ~ 0.38. They also found that this effect depends on 
the charge of the head group, and it disappeared when POPC 
(having a zwitterionic head group) was replaced by neutral 
dioleylglycerol (DOG). 

Our results in Fig. |3l suggest local extrema in {A ) when 
Xtap is between 0.16 and 0.5 (in addition to global minimum 
at Xtap — 0.5). In addition to the above-mentioned POPC 
and SS-1 study dSailv et all '2001'), similar and even more 
dramatic effects have been observe d for mixtures of POPC 
and sphingosine dSailv et all l2003l) . The existence of criti- 
cal concentrations in lipid membranes has also been theoreti- 
cally postulated by Somerharjuet al. ( Somerhariu et al., 1999; 
lyirtanenet al., 1998). In the present case, the local extrema 
for DMPC / DMTAP mixtures are within error bars, and there- 
fore may be interpreted as fluctuations of (A). To study such 
features in more detail one needs to decrease the fluctuations 
in the average area per lipid by, e.g., increasing the system 



size (^Lindahl and Edholml l200Cil) . This, however, is beyond 
the scope of the present study. 

Nevertheless, we decided to approach this issue from a dif- 
ferent perspective. To determine the average area per lipid 
separately for the two different components, we us ed the 
Voron oi tessellation technique in two dimensions ( Patr a et all 
l2003h . In Voronoi tessellation, we first calculated the center of 
mass (CM) positions for the lipids and projected them onto the 
x-y plane. A point in the plane is then considered to belong to 
a particular Voronoi cell, if it is closer to the projected CM of 
the lipid molecule associated with that cell than to any other 
CM position. As there is no unique definition for the area 
per molecule in a multi-component system, it is clear that the 
Voronoi results should be considered as suggestive rather than 
quantitative, providing insight mainly of the trends. 

Figure |3 demonstrates that the areas occupied by DMPC 
and DMTAP are distinctly different. For small xtap, the area 
per DMPC is considerably larger than that of DMTAP. For 
larger DMTAP mole fractions above xtap =0.5, the situa- 
tion is the opposite. This behavior is related to electrostatic 
effects and the ordering of acyl chains, and will be discussed 
in more detail in Sect. IIII.BI Here we only note that the fluc- 
tuations in {A) (see Fig.O at 0.1 < Xtap ^ 0.5 arise from 
fluctuations in the area occupied by DMTAP. Whether this is 
a true result due to, e.g., clustering of lipids in this region re- 
mains to be resolved. 



B. Ordering of lipid acyl chains 

Ordering of non-polar hydrocarbon chains in lipid bilay- 
ers is typically characterized by the deuterium order parame- 
ter ScD measured through •^H NMR experiments. If 6 is the 
angle between a CD-bond and the bilayer normal, the order 
parameter is defined as 

ScD^^icOS^e)-^ (1) 

separately for each hydrocarbon group. Since we employed a 
united-atom force field, the positions of the deuterium atoms 
are not directly available but have to be reconstructed from 
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FIG. 5 Deuterium order parameter IScol averaged over sn-1 and 
sn-2 chains for DMPC (top) and for DMTAP (bottom) as a function 
of DMTAP mole fraction. Small carbon atom numbers correspond 
to those close to the headgroup. 



the coordinates of three successive nonpolar hydrocarbons, 
assumi ng an ideal tetrahedral geometry of th e central CH2 
group jChiu et all fTool iHofsass et all l200 3f). In practice, 
we calcula ted Sg n following the standard approach described 
elsewhere frat ra et alll2003l) . 

Figure|5l shows l^cnl averaged over the two similar atoms 
in the sn-1 and sn-2 chains, for both DMPC (top) and DM- 
TAP (bottom) at different DMTAP concentrations. For the 
pure DMPC bilayer, we find IS'cdI ~ 0.18 close to the glyc- 
erol group of the molecule, in goo d agreement with recent 
experiments ("Petrache et al.', '200(f) and molecular dynam- 
ics simulation studies (Rog and Pasenkiewicz-Gierula, 20011 
ISmondvrev and Berkowitztl2001l) . 

The order parameter profile for the first seven hydrocarbons 
(from C2 to C8) is a kind of plateau, and the average value of 
ScD in this region, denoted by S'avc, is shown in Fig.|6] As 
expected, the results closely follow the change in the average 
area per lipid, i.e., a compression of the membrane is accom- 
panied by enhanced ordering of nonpolar acyl chains. Results 
of similar nature have been observed, e.g., in bilayer mix- 
tures of glycerophospholipids and cholesterol, where choles- 
terol both reduces the average area per molecule and enhances 
the ordering of lipid acyl chains (Hofsass et al., 2003). 

For pure DMPC the plateau value 5avc equals to 0.181 ± 
0.009, which is in very good agreeme nt with the exp erimen- 
tally measured value of 0.184 ( Petrac he et al I I2OOOI) . For a 
pure DMTAP bilayer, the plateau value of Scd is consider- 
ably smaller, being about 0.147 ± 0.011. Therefore chains in 
a pure DMTAP bilayer are on average more disordered than 
in a DMPC system, in agreement with our findings for {A ). 



C. Orientation of phosphatidylcholine headgroups 

Since the chemical structures of the acyl chains of DMPC 
and DMTAP are identical, it is obvious that the differences 




0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 



FIG. 6 Plateau order parameter Save calculated by averaging Scd 
over C2 to C8 hydrocarbons. Shown are Save for DMPC (solid lines 
with solid circles) and DMTAP (dashed lines with open circles). 
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FIG. 7 Results for the probability distribution function P(a) vs the 
angle a between the P-N vector (of DMPC headgroups) and the 
bilayer normal. 



between their behavior are due to their headgroups. Since the 
headgroup of DMPC is zwitterionic (cf. Fig.Q, it possesses 
a dipole moment along the P-N vector. The electrostatic po- 
tential across a monolayer thus depends sensitively on the dis- 
tribution of the angle a between the P-N vector and the inter- 
facial normal n (where n has been chosen to point away from 
the bilayer center along the z coordinate). 

Figure shows the probability distribution function P{a) 
for the angle in question. For a pure DMPC bilayer we find 
the distribution to be wide, thus allowing the P-N vector 
to fluctuate substantially, pointing at times in the direction 
of the membrane normal as well as into the bilayer interior 
The average angle found in this case is about (80 ±1)° (see 
Fig. |8}, i. e., the PC heads are on average almost parallel to 
the membrane surf ace. This is in agreement with experimen- 
tal o bservations (iHauser et all fl981: Scherer and Seelij, 
Il989l) as well as with recent computer simula- 
tions ("Gabdoulline et al.', '199(?, 'Pasenk iewicz-Gierula et all 
1999; Smondvrev and Berkowitz, 1999). 

Figure further reveals the major role played by DMTAP. 
The addition of even a small amount of DMTAP leads to a pro- 
nounced re-orientation of the P-N dipole vector. As xtap is 
increased, the profile of the distribution becomes considerably 
narrower and its maximum shifts to smaller angles. This trend 
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FIG. 8 The average angle (a) between the P-N vector of DMPC 
and the bilayer normal, shown as a function of DMTAP mole frac- 
tion. Full circles: Results averaged over all DMPC molecules in a 
given system. Open circles: Results averaged over only those DMPC 
molecules that are beside DMTAP. See text for details. 



continues up to the high-concentration limit xtap ~ 0.75, be- 
yond which the distribution is essentially similar with the case 
found for xtap — 0.75. 

Results for the average angle (a) between the P-N vector 
and the membrane normal shown in Fig.|8]are consistent with 
this picture. On average, upon increasing xtap, PC head- 
groups become more and more vertically oriented. Also this 
has been observed in se veral experiments ( Scherer and Seelii, 
fT989HSailv et al.Ll2003li2()0lUZantl et all ll999b>. Moreover, 
our findings are in fairly good agreement with an atomistic 
MD study of a compl ex comprised of DNA and a mixture of 
DMPC and DMTAP (B andvopadhvav et all \l99% . in which 
the average angle between the P-N dipole vector and the bi- 
layer normal was found to be (50±8)° at an almost equimolar 
mixture of DMPC and DMTAP. In the present case without 
DNA, we found (a) (42 ± 2)° at xtap = 0.5. 

Perhaps surprisingly, the correlation between the average 
area per lipid (Fig.|3} and the re-orientation of the P-N dipole 
is not complete. As Fig. |8l shows, the re-orientation extends 
by and large linearly up to xtap — 0.75, while the membrane 
compression completes at y t\v ^ 0.5. This contradicts the 
conclusions of Saily et al. dSailv et all I2OOI ) who studied 
POPC/SS-1 cationic lipid mixtures using the Langmuir bal- 
ance technique and suggested that the maximal average an- 
gle between the P-N dipole vector and membrane surface is 
achieved at the cationic lipid concentration that corresponds 
to the point where the membrane compression ends. 

A closer inspection of Figs.|3and|S]shows that the observed 
reduction of the average area per molecule is likely related to 
the re-orientation of PC dipoles, and this in turn is related to 
the role of electrostatic interactions between DMPC and DM- 
TAP headgroups. To bridge the two issues, we propose the 
following schematic scenario. At small xtap where the DM- 
TAP molecules are far apart and their mutual interaction is 
rather weak, we essentially suggest that the role of DMTAP is 
to re-orient the headgroups of those DMPC molecules that are 
beside a DMTAP molecule. This favors more dense packing at 
small Xtap, leading to a reduction in {A), and consequently 



to a minimum in the area per molecule at intermediate concen- 
trations as for large xtap the repulsive electrostatic interac- 
tions between TAP headgroups enforce {A ) to be expanded. 

To validate this scenario, we complemented our results 
in Fig. by calculating the probability distribution function 
P{a) for those DMPC molecules that are nearest neighbors 
to DMTAP 

As a criterion that a DMPC and a DMTAP form a pair, 
we monitored the distance between PC phosphorus and TAP 
nitrogen. For that, we first calculated the radial distribution 
functions (RDFs) between pairs of Ppc and Ntap and deter- 
mined the distance Tnn at which the RDF had its first mini- 
mum after the main peak (see also Sect. IIII.Fl . The distance 
obtained in this fashion (r„„ « 0.665 nm) (and found not to 
depend on xtap) was applied to identify the DMPCs residing 
next to a DMTAP. As shown in Fig.|8] the re-orientation of the 
P-N vector of these DMPCs is considerably stronger at small 
DMTAP concentrations as compared to that averaged over all 
DMPC Hpids. 

The above results imply that at small xtap the effect of 
DMTAP on the re-orientation is mainly local, i.e., the alter- 
nating PC and TAP headgroups pack more tightly than in a 
pure DMPC system. This idea is supported by the results 
for the radial distribution functions discussed in Sect. IIII.FI 
Beyond the small-xxAP regime, for intermediate concentra- 
tions 0.3 ^ Xtap 0.5, further increase in the concentration 
of DMTAP continues to increase the number of units com- 
posed of PC and TAP heads, thus favoring a reduction in {A). 
However, as repulsive electrostatic interactions between DM- 
TAP molecules also become more and more important, the 
two effects compensate each other and {A ) is found to be ap- 
proximately constant. Finally, for large xtap, the repulsive 
electrostatic interactions between TAP groups dictate the case 
discussed here and lead to an enhancement of the average area 
per molecule. 

Though this picture does not account for the explicit influ- 
ence of counter-ions, it grasps essence of the process. The 
effect of counter-ions is discussed separately in Sect. lIII.Fl 



D. Density profiles of lipid headgroups and chloride ions 

To quantify the locations of charge groups and counter-ions, 
we computed the density profiles across the bilayer, separated 
into the different constituents of the system. The positions 
of all atoms in the system were determined with respect to the 
instantaneous center of mass position of the bilayer, exploiting 
mirror symmetry such that atoms with z < were folded to 
z > (the bilayer center being at z = 0). 

Figure |9] shows the scaled number densities pn{z) for a 
few selected cases. Additionally, we note that the essential 
information is given by the positions of the density maxima 
depicted in Fig. 110! 

At small Xtap, the density maxima of the nitrogen and 
phosphorus atoms in the DMPC heads almost coincide (see 
Fig. |9}. The density profile of nitrogen in DMPC is never- 
theless broader and extends further out of the bilayer plane. 
For larger xtap, the density profiles of phosphorus and ni- 
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z. [nm] 

FIG. 9 Scaled number densities pn{z) for three DMPC/DMTAP 
mixtures with xtap = 0.0 (top), xtap = 0.5 (middle), and 
Xtap = 1.0 (bottom). The case z = corresponds to the center 
of the bilayer. 



trogen are distinctly separated, and nitrogen in particular ex- 
tends rather deeply into the water phase. The TAP group rep- 
resented by the nitrogen atom, however, is found to be deep in 
the bilayer It seems obvious that these two issues are related, 
i. e., the density profiles of nitrogens in PC and TAP groups 
are well separated due to the electrostatic repulsion that es- 
sentially leads to the re-orientation of PC headgroups. These 
results are hence consistent with those in Sect. lIII.Cl and reflect 
the dependence of DMPC headgroup orientation on xtap ■ 

Interestingly, while being attracted by the DMTAP head- 
groups, the chloride anions cannot penetrate the outer bound- 
ary of the bilayer formed by the DMPC choline groups. This 
is in a sense to be expected since the DMPC headgroup is 
longer than the TAP group and thus extends further outward 
from the bilayer. There is thus a significant amount of shield- 
ing of the chloride ions in the presence of DMPC. Only for an 
almost pure DMTAP layer are the chloride ions located in the 
vicinity of the TAP headgroups. 



E. Charge density, electrostatic potential, and orientation 
of water 

The charge distribution shown in Fig.^Jwas calculated in 
the same fashion as the density profiles. The results are clearly 
reminiscent of the density profiles in Fig. |9] and demonstrate 




z [nm] 

max 

FIG. 10 Maxima of the density profiles, Zmax, for phosphorus and 
nitrogen atoms from the DMPC headgroups, and of the density pro- 
file of chloride ions. The maxima are shown for the z coordinate in 
the direction of the membrane normal, shown as a function of DM- 
TAP molarity. The dashed line marks the position of the membrane- 
water interface determined from the condition that the densities of 
water and lipids (in Fig.|9) are equal. 



the competition between charged PC and TAP groups, CI an- 
ions, and water The role of the TAP group is prominent. 

From the charge densities, we also computed the electro- 
static potential across the bilayer. The results are shown 
in Fig. where the potential at the center of the bi- 
layer has been set to zero. For a pure DMPC bilayer 
we obtain —0.578 V, which is in agreemen t with previous 
MD simulation studies fchiu et al.', 1995: .GabdouUine et all 
1 19961 FSmondvrev and Berkowitz. 1999). Experimental data 
for phospholipid membranes ranges from — 0.2V to —0.6V 
tiClai-ke. 2001: Flewelfing and HubbeL 1986:«Gawrisch et all 
II992l:lHladkv and HavdonIll973HPickard and BenzLll978h . 

The charge of the DMTAP headgroup is mainly compen- 
sated by, depending on the value of xtap, chloride ions or 
DMPC phosphate groups. Most of the electrostatic potential 
across the bilayer thus is not due to the DMTAP itself but 
rather due to the re-orientation of the DMPC headgroups. A 
clear indication of this is that the potential build-up saturates 
at Xtap — 0.75, i. e., at the same value at which the distri- 
bution of headgroup orientation saturates (see Fig.|8}. The to- 
tal potential of the bilayer increases with increasing DMTAP 
concentration, with a difference of 0.6 V between pure DMPC 
and pure DMTAP. This increase agrees well with the experi- 
menta l data on cationic POPC / SS-1 monolayers ( Sailv et all 
I2OOII) . and the authors offer the same explanation for their ob- 
servations. 

Many of the conclusions drawn from the charge den- 
sity already follow from the number densities presented in 
Sect. IIII.DI since for charged particles number density and 
charge density are trivially related. Water, however, has an 
additional internal degree of freedom, and a quick discus- 
sion of the orientation of the water molecules seems appropri- 
ate. As seen from Fig.[Oj the average direction of the water 
dipoles in the membrane-water interface region is inverted for 
Xtap = 0.5 Xtap = 1-0. This is closely related with the 
familiar "hump" close to the interface, which is due to a sub- 
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FIG. 1 1 Charge densities p{z) across a single leaflet for xtap equal 
to 0.0 (top), 0.5 (middle), and 1.0 (bottom). The case 2 = corre- 
sponds to the center of the bilayer. Charge densities are shown as 
solid lines. In addition, the component-wise contributions due to 
DMPC (full circle), DMTAP (open circle), water (dashed line), and 
chloride ions (star) are displayed. To reduce the noise i n the data, the 
charge densities shown here were first fitted to splines jThiisse et all 
[l998l) . The error bars are of the same size as the symbols. 




FIG. 12 Electrostatic potential V{z) across cationic bilayers at dif- 
ferent DMTAP mole fractions. 



tie imbalance between the orientation of the water molecules 
and hpid headgroups dChiu et allll995l) . At higher DMTAP 
concentrations this "hump" disappears. A related issue con- 
cerns the pure DMTAP bilayer, in which case the density pro- 
file of water penetrates rather deep into the membrane (see 
Fig.|9}, extending up to the interface region between the polar 
TAP group and the hydrophobic core. This is in accord with 
the interpretation of Fourier-tr ansform infrared sp ectroscopic 
measurements by Lewis et al. dLewis et al l l200li) . 




2.0 2.5 

z [nm] 

FIG. 13 Projection of the water dipole unit vector il{z) onto the 
interfacial normal n, yielding P{z) = {f^{z) ■ ft) = {costf)}. Here 
z — corresponds to the center of the bilayer and the bilayer normal 
n is chosen to point away from the bilayer center. 



F. Radial distribution functions and coordination 
numbers 

To characterize the structure of the membrane-water inter- 
face region in more detail, we computed various radial distri- 
bution functions (RDFs) between the atoms in the headgroups, 
CI, and the oxygens. Here we briefly discuss the most relevant 
results. 

The RDFs between the center of mass positions indicated 
that the leading (main) peak for DMPC -DMPC and DM- 
TAP -DMTAP pairs was rather broad and at about 1.0 nm 
(data not shown). For DMPC -DMTAP pairs, however, the 
main peak of the RDF was much closer, at around 0.7 nm. 
This supports the conclusion made in Sect. lIII.Cl i.e., DMPC 
and DMTAP form units that allow more dense packing than 
in a pure DMPC bilayer. 

The Npc-Npc and Ntap-Ntap pairs were found to 
be rather far apart, the position of their main peak being at 
about 0.83 nm, while the Npc -Ntap pair was slightly closer 
(0.8 nm). The positions of the main peaks did not depend on 
Xtap ■ As for the RDFs of the phosphorus atoms in the PC 
headgroups, its main peak with respect to Npc and Ntap was 
found to be at a much closer distance, at 0.465 nm for Npc 
and 0.485 nm for Ntap- Again, the positions of these peaks 
did not depend on the DMTAP concentration. 

We also calculated the coordination numbers for the phos- 
phorus and nitrogen atoms at different DMTAP concentra- 
tions. These are shown in Fig. 1141 (top). It turns out that in 
the range from xtap = to 0.5 the PC nitrogens are to an 
increasing extent being replaced by Ntap in the vicinity of P. 
This has twofold consequences: First, the electrostatic attrac- 
tion between N+ (TAP) and P^ (PC) enhances the compres- 
sion of the bilayer for 0.0 < Xtap ^ 0.5 (see Fig.|3jl. Second, 
the decreasing coordination number for P-Npc with xtap 
support the view that the DMPC nitrogens are pushed towards 
water, thereby PC headgroups are re-oriented to a more verti- 
cal alignment with respect to the membrane plane (Fig.|8}. 

We conclude this section with a discussion of the location 
of chloride counter-ions. The positions of the main peaks of 
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FIG. 14 Top: Coordination numbers Nc of DMPC phosphorus with 
DMPC nitrogen (solid line with solid circles) and with DMTAP ni- 
trogen (dashed line with open circles) plotted versus xtap- Bottom: 
Coordination numbers Nc of DMPC nitrogen (solid line with solid 
circles) and DMTAP nitrogen (dashed line with open circles) with CI 



the RDFs for both types of nitrogens (in PC and TAP groups) 
with respect to CI ions are identical, about 0.475 nm, and do 
not depend on xtap (data not shown). In Fig.[2](bottom) we 
plot the coordination numbers of chlorides in the vicinity of 
both types of nitrogens as a function of DMTAP concentra- 
tion. The figure confirms that CI ions are preferentially bound 
to PC nitrogens rather than to Ntap- This holds up to a DM- 
TAP mole fraction of about 0.75. The explanation for this is 
straightforward: As xtap increases, the PC headgroups be- 
come more and more vertically oriented with respect to the 
bilayer plane. This, in turn, makes PC nitrogens more easily 
accessible for the CI ions. In contrast, small TAP heads are 
located much deeper in the membrane surface region than the 
PC heads and therefore are able to attract fewer chlorides re- 
gardless of the fact that TAP heads carry a net positive charge. 
Interestingly, when the re-orientation of PC heads is accom- 
plished (at Xtap ~ 0.75) the coordination number for Npc - 
CI pairs seems to saturate, see Fig.ll4l(bottom). 



IV. SUMMARY AND CONCLUSIONS 

Drug delivery and gene therapy have attracted substantial 
interest due to their importance in treating human diseases. 



As far as experimental work is concerned, particular attention 
has been paid to non- viral delivery vectors such as cationic li- 
posomes characterized by a number of desired properties such 
as high efficiency and lack of toxicity. Consequently, it is sur- 
prising how little is known about the atomic-level details of 
cationic lipid bilayers. Essentially, this is due to the lack of 
molecular simulations of these systems, the study by Bandy- 
opadhyay et al. ( Bandvopadhvav et al., 1999) being the only 
exception, to the authors' knowledge, in this regard. 

As a first step toward a detailed understanding of cationic 
membrane-DNA complexes on atomic level, we have em- 
ployed extensive molecular dynamics simulations of lipid bi- 
layer mixtures composed of cationic DMTAP and neutral 
(zwitterionic) DMPC. Such binary DMPC / DMTAP mixtures 
have been studied widely through experiments, and have been 
shown to form stable complexes with DNA (lArtzner et all 
ll998l:lPohIe et alll2000HZantI et allll999albh . In the present 
work, we have focused on the influence of the composition of 
the cationic bilayer on its structural and electrostatic proper- 
ties. For this purpose we studied numerous DMPC /DMTAP 
mixtures in the liquid-crystalline phase by varying the mole 
fraction of DMTAP, xtap, from the pure DMPC to the pure 
DMTAP bilayer 

We have found that the properties of the DMPC /DMTAP 
bilayer mixture are largely dominated by the electrostatic 
properties of the headgroup region around the membrane- 
water interface. Most notably, our results indicate that there 
is a strong interplay between the PC and TAP groups together 
with the CI counter-ions that concentrate in the vicinity of the 
bilayer-water interface. 

The interplay between the PC and TAP groups leads to 
a number of intriguing observations. The key factor here 
is the re-orientation of PC groups due to an introduction of 
DMTAP in the bilayer. The re-orientation of the PC head- 
groups arises from electrostatic interactions that lead phos- 
phate and choline groups to rearrange their positions with re- 
spect to the cationic TAP. This effect is enhanced as xtap is 
increased, and extends up to large molar fractions of approxi- 
mately Xtap = 0.75. Beyond this limit a further increase of 
DMTAP concentration has no additional effect on the orienta- 
tion of PC headgroups. Interestingly, at small xtap the effect 
of the re-orientation is of local nature, i.e., the P-N dipoles 
of DMPCs beside DMTAP molecules re-orient considerably. 

At small molar fractions of DMTAP, the re-orientation of 
PC dipoles leads to considerable compression of the bilayer as 
alternating PC and TAP groups are able to pack more tightly 
than in a pure DMPC bilayer. The minimum of the area per 
lipid at Xtap ~ 0.5 is about 12 % smaller than in the pure 
DMPC bilayer. A further increase of xtap leads to major 
expansion of the bilayer This is essentially due to an in- 
creasing number of neighboring TAP groups whose cationic 
nature leads to repulsive electrostatic interactions that do not 
favor close packing. As expected, the ordering of acyl chains 
closely follows the change in the area per lipid. When these 
results are summarized, the present view of the average area 
per lipid coupled to the re-orientation of the headgroups can 
be summarized schematically as in Fig. 1151 

In view of future studies of DNA-membrane systems, it 
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FIG. 15 A proposed schematic picture of the observed change in 
the area per lipid vs xtap. Only headgroups of the lipids are shown 
here, and water (not shown) is above the headgroups. 



is important to pay attention to the influence of DMTAP on 
the electrostatic properties of the membrane, including the 
increase in the electrostatic potential across the bilayer and 
the ordering of water in the vicinity of the membrane-water 
interface. Another often ignored aspect of electrostatics is 
that th e ionic buffer hquids m ay affect membranes signifi- 
cantly jBockmann et alll2003l) . 

Perhaps the most significant observation in this study is the 
spatial re-arrangement of PC and TAP headgroups which is 
expected to play a significant role in the condensation of DNA 
onto the membrane surface. The cationic TAP and choline 
groups then play a key role as the anionic phosphate groups 
of DNA come into contact with the membrane. While the 
present study clarifies some of the underlying questions re- 
lated to binary mixtures of cationic and neutral (zwitterionic) 
lipid membranes, further atomic-level studies are essential to 
resolve other important issues related to DNA-membrane sys- 
tems, such as the influence of salt and its screening effects, 
and the stability and interface properties under those condi- 
tions. Work in this direction is under way. 
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